school_plot <-
  read.csv("./_3_data/indonesia_school_data.csv") %>%
  mutate_at(vars(num_high_school, num_priv_high, num_mid_school,num_priv_mid, num_prim_school,num_priv_prim), funs( . = ./(pop/1000))) %>%
  select(pop, num_high_school_., num_priv_high_., num_mid_school_., num_priv_mid_.,  num_prim_school_.,num_priv_prim_.) %>%
  pivot_longer(cols = num_high_school_.:num_prim_school_., names_to = "type", values_to = "value") %>%
  mutate(type_label = case_when(str_detect(type, "priv") ~ "Private",
                                TRUE ~ "Public"),
         level = case_when(str_detect(type, "high") ~ "High School",
                           str_detect(type, "mid") ~ "Middle School",
                           str_detect(type, "prim") ~ "Primary School")) %>%
  filter(type_label == "Public") %>%
  ggplot(aes(x=log10(pop), y = value)) +
  geom_point(shape = 21, alpha = 0.3, fill = "lightgrey") +
  geom_smooth(color = "black") +
  facet_wrap(level ~ ., ncol = 3, scales = "free_y") +
  stat_summary_bin(fun='mean', bins=20, fill = "lightgrey",
                   shape = 21, size=2.5, geom='point') +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), 
        panel.grid.major.x = element_blank(),
        axis.line.y.left = element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(),
        legend.title = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank()) +
  ylab("Schools per 1000 Residents") +
  xlab("Population (logged)")

ggsave("./_4_outputs/figure_oa_i.pdf", plot = school_plot, width = 10, height = 4)
